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Abstract 



We present a new combinatorial method for the calculation of the nuclear level density. 
It is based on a Monte Carlo technique, in order to avoid a direct counting procedure which 
is generally impracticable for high-A nuclei. The Monte Carlo simulation, making use of 
the Metropolis sampling scheme, allows a computationally fast estimate of the level density 
for many fermion systems in large shell model spaces. Wc emphasize the advantages of 
this Monte Carlo approach, particularly concerning the prediction of the spin and parity 
distributions of the excited states, and compare our results with those derived from a 
traditional combinatorial or a statistical method. Such a Monte Carlo technique seems 
very promising to determine accurate level densities in a large energy range for nuclear 
reaction calculations. 
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1 Introduction 



The density of nuclear levels provides information about the structure of highly excited nuclei, 
and is also a basic quantity in nuclear reaction theory. For many years, measurements of 
nuclear level densities have been interpreted in the framework of an infinite Fermi gas model 
(see e.g. Jll), and most calculations have employed the well-known Bethe formula since it 
predicts a general trend in reasonable accordance with experimental data (i.e., an approximately 
exponential increase of the level density with both the excitation energy and the mass number). 
The Bethe formula is based on the statistical model (see e.g. |0, ^), and involves three major 
assumptions: (z) the independent particles assumption, allowing to write a nuclear partition 
function in a simple form using single-particle energies, (ii) an equidistant spacing of the single- 
particle states near the Fermi level, and (iii) the saddle-point approximation used to calculate 
the inverse Laplace transform of the partition function. As a result of assumptions (i) and (ii), 
this simple model includes neither odd-even effects nor shell effects, and it has been pointed 
out in ref. P| that assumption {in) is not very good in the nuclear case. Therefore, various 
phenomenological modifications of this simple formula have been proposed for use in practical 
calculations to account for these effects (see ref. for a recent review). These are mostly 
based on the ad-hoc hypothesis that the same functional form of the energy dependence as in 
the equidistant-level case is valid, some free parameter(s) being introduced in order to match 
experimental results. In the Newton- Cameron shifted Fermi-gas model 0, |^, the odd-even 
effects are included by means of a pairing energy shift. The effective excitation energy is 
reduced by the conventional pairing energy for the odd-mass and even-even nuclei, resulting 
in a lower level density. The pairing energy is determined from experimental odd-even mass 
differences, and the level density parameter a is the only adjustable parameter. This model is 
able to describe experiments in a narrow energy interval around the neutron binding energy, 
most experimental data being extracted from the neutron resonances, but it does not allow to 
extrapolate the value of the level density in either the low- or the high-energy regions. This 
means that shell effects are not properly taken into account. 

In order to solve this problem, Gilbert and Cameron |^ proposed another formula (the com- 
posed four-parameter formula) which combines the shifted Fermi gas formula at high excitation 
energies with a constant temperature formula (see e.g. |TD[) for lower energies. By fitting the 
four constants in both regions, experimental data may be well reproduced. Another simple ap- 
proach, the back-shifted Fermi gas model was later proposed ||Tl|, |13|, [I^ in order to account 
for shell effects. In this model, both the energy shift and the level density parameter a are con- 
sidered as adjustable parameters, which allows to obtain a reasonable fit to the experimental 
level densities over a wider range of excitation energies. Afterwards, some phenomenological 
methods have been proposed to correctly describe the thermal damping of shell effects with 
increasing excitation energy (e.g. |TB|, |T^, 0). In those methods, the idea is to reproduce as 
well as possible the energy dependence of the level density parameter a taken from microscopic 
calculations, in order to give better absolute values of the nuclear level density. 

However, most of the above-mentioned semiempirical approaches are based on various dras- 
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tic approximations, and their shortcomings at matching experiments are often overcome only 
by parameter adjustments. Therefore, the quahty of the phenomenological prescriptions for a 
wide range of mass numbers can not be taken as a support for their accuracy outside the narrow 
regions of excitation energy (~ 5-15 MeV) and angular momentum (~ 0-5 h) to which most 
of the experimental knowledge of level densities is confined. Numerical statistical calculations 
have also been developed to calculate a more realistic level density using the single-particle level 
scheme of the shell model P, and the BCS formalism has been included ||T9| , pO[| in order 
to account more properly for pairing effects. However, those microscopic calculations are still 
depending on the assumptions of the statistical model, namely the saddle-point approximation 
which is not very good at low energy. Moreover, the spin and parity distributions raise unsolved 
questions and call for improvement. 

The advent of high speed computers has made possible to use methods of calculation which 
do not depend on closed formulae, such as the combinatorial method [^1], 0, The level 



density is then calculated numerically by performing an exhaustive counting of the nuclear 
excited configurations. It yields an exact - at least in the independent-particle picture of the 
nucleus - level density, but is very time consuming and becomes intractable at high excitation 
energies or for large shell model spaces (i.e. for high-A nuclei). The pairing interaction may 
be included in the calculations by applying the BCS theory, but at the expense of a strongly 
increased computation time. More recently, combinatorial calculations have been carried out 
to estimate level densities with fixed exciton numbers |2^, for use in preequilibrium models. 



A second method of computing nuclear level densities directly from a set of single-particle 



states has been proposed by Williams |2^, in order to avoid the direct enumeration and clas- 
sification of states which renders the combinatorial method generally cumbersome. It is based 
on the repetitive use of recursion relations to expand the grand partition function, but it does 
not rely on the saddle-point approximation like the traditional statistical method. This so- 
called recursive method provides an exact calculation of level densities, but, compared with 
the combinatorial approach, it has the advantage of requiring far shorter computation time. 



It has been been extended in ref. [^] in order to describe the spin-parity distribution, and 
also to provide state densities with fixed particle-hole number. Unfortunately, it cannot treat 
residual interactions exactly, at least in its classical form, so that it is useful in the framework 
of a non-interacting Fermi gas model only. Spectral distribution methods ^ (based on 
moment techniques to calculate averaged distributions of the Hamiltonian eigenvalues) have 
been used to provide level densities taking residual interactions into consideration. However, 
these methods are generally limited to a low-order expansion, which is not very accurate for 
large shell model spaces. Thus, for different reasons, both the combinatorial and the recursive 
methods are inappropriate to yield microscopic level densities including residual interactions in 
a reasonable computation time. 

In this paper, we propose a Monte Carlo method as an alternative to the traditional combi- 
natorial method. We calculate the level densities of spherical nuclei by means of a Monte Carlo 
simulation based on the Metropolis sampling scheme |2^. We show that this technique, initially 
devised for the simulation of molecular systems, can be successfully applied to the nuclear level 
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density problem. It yields an exact - save the statistical errors - level density within a short 
cpu time, whatever the mass number of the nucleus or the considered excitation energy. We 
also use this Monte Carlo method to calculate the spin and parity distributions of the excited 
levels. In this first attempt, we restrict ourselves to the case of spherical nuclei, but the method 
could be easily extended to deformed nuclei. A detailed description of the Monte Carlo method 
is given in Section |[ In Section ^, we check the reliability of this Monte Carlo procedure by 
comparing our results with those of a direct counting and with microscopic statistical calcu- 
lations. Section ^ is devoted to the inclusion of pairing interaction in this model. Finally, we 
conclude in Section |^. The Monte Carlo implementation of the spin coupling and the statistical 
errors originating from the method are discussed in the Appendices. More detailed results and 
a comparison with experimental data are presented in a second paper [HO . 



2 The Monte Carlo method 



As noted first by van Lier and Uhlenbeck , the nuclear level density problem in its simplest 
form (i.e. the uniform spacing model with one kind of particles) is equivalent to a well-known 
combinatorial problem in number theory, the partitioning of integers, first solved by Euler in 
the eighteenth century. Indeed, for a system of N fermions, the number M = gU represents 
the number of quanta which are to be distributed among the N fermions, g and U being the 
single-particle state density and the excitation energy, respectively. Thus, the state density 



u;{U)=gpN{M) 



(1) 



is proportional to the number pn{M) of partitions of the integer M into no more than 

p^] ). Note that an approximate solution for Euler's problem was derived 



positive parts (see 
by Hardy and Ramanujan in ref. 
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Pn{M) 



exp(7r 



|M) 



(2) 



when N > M (i.e., when the excitation energy is not sufficient to excite a particle from the 
lowest orbit, which is clearly the case for an infinite model). This approximation is in fact 
equivalent to Bethe's formula 

exp(2-\/af7) 



(3) 



with the parameter a = ^g. As investigated by Euler, the exact pn{M) can be calculated by 
means of the recursive relation 



Pn{M)=pm-i{M)+pm{M-N) . (4) 
It is instructive to notice that this relation is of the same class as Williams recursive relation 



25 1 for the calculation of level densities in the uniform spacing model. However, the Williams 
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method is also applicable for a non-equidistant single-particle level scheme by expressing the 
single-particle energies as multiple of an energy unit (chosen approximately equal to the accu- 
racy on the energies, e.g. 0.01 MeV). In this case, the level density problem comes to a more 
complex problem known in computer science as the subset-sum problem, for which the best 
algorithmic solution is a recursive one [Q. Unfortunately, as mentioned above, the recursive 



method can not incorporate residual interactions, so that recourse to a direct counting method 
cannot be avoided when treating pairing interaction. On the other hand, it is well known that 
a direct counting of all the excited nuclear levels is a very tedious task, almost impractica- 
ble on the currently available computers whenever one is interested in medium mass nuclei at 
intermediate energies. Moreover, a direct counting procedure cannot easily treat the residual 
interaction in nuclei because it is then exceedingly time-consuming. Thus it is natural to resort 
to a Monte Carlo technique in order to avoid this exhaustive counting of the excited levels. 
Furthermore, it is well known that Monte Carlo methods provide generally very efficient algo- 
rithms for solving combinatorial problems (e.g. combinatorial optimisation problems solved by 



simulated annealing |]35|, |3^). This is the main reason why we investigate the use of a Monte 
Carlo technique in the context of the nuclear level density problem, and this idea has been first 
explored in ref. This Monte Carlo method becomes essential if the residual interactions 



are taken into account, since the recursive method is intractable in this case. For simplicity, 
however, we consider the nucleus as a system of non-interacting fermions in this Section. This 
choice will allow a comparison with a direct counting in Section ^ The inclusion of the residual 
interactions is the subject of Section ^. 

2.1 Principle of the Monte Carlo procedure 

We present in this Section a Monte Carlo method based on the Metropolis sampling scheme 
which makes possible the determination of a combinatorial level density without direct counting. 
The Metropolis technique was devised for the simulation of molecular systems p9|; it has since 



been widely applied in statistical mechanics (e.g. to calculate thermal averages) |3^, ^ ^ 
We show here that it is convenient for describing nuclear level densities as well. Note that 
Monte Carlo simulations related to the level density problem have been performed in ref. E 



to obtain quantities of interest for the determination of the parity distribution. However, that 
work was done in the framework of the statistical model, which appears to be unable to treat 
correctly the parity dependence of the level density ^ . 



The Monte Carlo simulation is based on a random sampling of a very small fraction of 
the excited states in the considered range of excitation energy ||4^. The resulting sample is 



assumed to be representative of the whole configuration space, in analogy with what is done 
when estimating a multi-dimensional integral by a Monte Carlo procedure (see e.g. [^). Since 
the state density of a system of non-interacting fermions may be trivially expressed as 

uj{U, J,7t) = J2 SiU - Uc) 6j,j^ 5.,^^ , (5) 
c 
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i.e. as a mult i- dimensional summation in a discrete space of configurations (each configuration 
C having an excitation energy Uc, spin Jc, and parity ttc), it can be efficiently obtained by a 
Monte Carlo sampling procedure. Note that each iV-particle state is called here a configuration; 
it corresponds to the set of occupation numbers for all the single-particle states. (This is 
often called a combination in combinatorial methods terminology). It is clear that such a 
procedure cannot yield the whole ensemble of 5-spikes, but it is not important in practice, 
since only the average state density in energy bins is needed. Thus, the properties of the whole 
spectrum of excited states (i.e., energy, spin, parity, or other quantities) can be simply derived 
by extrapolating from the sample and applying the appropriate scale factor. This procedure 
provides an exact state density in the sense that it is independent of the grand-canonical 
statistical formalism (i.e., the energy and the number of neutrons or protons are not fixed only 
on average, which is a well-known source of errors). Of course, one has to keep in mind that 
the results are inherently affected by a statistical error. 

A Monte Carlo estimator for the state density can be simply written as 

uj{U, J, tt) ~ 5* a;smp(f^, J, tt) , (6) 

where S is the scale factor (to be determined) and cjsmp stands for the state density in the 
Monte Carlo sample. One has 

N 

^^smp{U, J,n) = ^6{U - Ui) 6j^j^ 5^,^, , (7) 

i=l 

where the Cj are randomly (uniformly) sampled configurations (of energy Ui, spin Jj and parity 
TTj) in the whole space, and N is the number of configurations in the Monte Carlo sample. The 
scale factor can be expressed as 

5 = (E 1)/^ > (8) 
c 

but it can not be calculated by this method, since the summation is prohibitively long in our 
configuration space {J2c 1 is intractable). It has to be determined by an independent method 
(see below). 

Such an approach to the nuclear level density problem, introducing random excited config- 
urations, may be quite justified. Indeed, the probability theory has already been applied with 
success in the statistical model to determine the asymptotic gaussian-shaped spin distribution 
and the parity distribution (see e.g. [§). Furthermore, the nuclear level density problem may 
be seen as the problem of exploring a very large multidimensional space of configurations. A 
Monte Carlo simulation is recognized to be the best technique for solving this problem, and 
is in fact the only available one in many cases Note that, although the method is exact 
in principle, the resulting cu is inherently affected by a statistical error which scales like 1/\/N 
(see Appendix ^ . The counterpart of this problem is that the accuracy of the results can be 
imposed by choosing an appropriate size N for the sample, and do not depend on the actual 
number of states in the considered energy range. (For instance, a rough estimate may be ob- 
tained within a very short computation time by sampling only a few configurations.) On the 
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contrary, when doing a direct counting of the states, the whole configuration space has to be 
explored, yielding an exact result within a computation time which can be prohibitively long. 
There is no other alternative. This fact constitutes the major advantage of a Monte Carlo 
procedure. 

2.2 The Metropolis sampling scheme 

It is clear that a simple random sampling of the states would be completely unefficient, since 
the state density uj{U, J, n) increases exponentially with excitation energy U. It is necessary 
in our simulation to hinder the random sampling of the highly excited configurations which 
are by far the most numerous, in order to achieve the same level of accuracy over the whole 
considered energy range. Therefore, we have to apply an importance sampling method (see e.g. 
[p8|). Let us introduce a (discrete) probability distribution P{C). The Monte Carlo evaluation 
of u;(?7, J, tt) becomes 



whenever each Ci is randomly sampled with a probability P{Ci). The problem is that the 
probabilities P{Ci) cannot be simulated (neither calculated) because it is impossible to perform 
the normalization in our configuration space. The Metropolis method provides a solution for 
this problem by only referring to ratios of probabilities, and makes thus an importance sampling 
procedure practicable. 

The Metropolis sampling scheme allows to produce any random variable in many dimensions 
with a given probability distribution. It only requires the ability to calculate a weight function 
(proportional to the probability) for a given value of the integration variable. Let us define 
a weight function W{C) that represents in our problem the relative (unnormalized) sampling 
probability of the configuration C, i.e., W{C) oc P{C). At this point, W{C) is arbitrary, 
but we will show how to define it properly in the following. The Metropolis method is based 
on a guided random walk which proceeds through configuration space according to a given 
matrix of transition probabilities p. Each element p{Ca — > Cb) of this matrix corresponds to 
the probability of performing a step from the state (or configuration) a to the state h. The 
Metropolis method prescribes the choice of this stochastic matrix p in order to obtain a random 
walk with a limit distribution (i.e., the distribution of the frequencies with which each state 
occurs when letting the random walk run to infinity) equal to a given P{C): 



where q is the so-called underlying matrix of the random walk, and will be defined below. In 
order to completely define the random walk process, the diagonal elements of the stochastic 






(10) 
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matrix must be taken as 

piCa ^Ca) = l- J2PiCa ^ C,) . (11) 

This Metropolis prescription has been proved to be rigorous for a discrete space which is 



just the case for our nuclear configuration space. Indeed, it is easy to check, using equation (jlQ]), 
that the detailed balance P{Ca) p{Ca Ch) = P{Cb) p{Cb Co) is verified provided that q is 
a symmetric matrix. This clearly implies that, conditionally on the ergodicity of the random 
walk, the correct limit distribution is reached. As already mentioned, expression (|ToD only 
refers to the weight function W{C), and is thus tractable. 

From the definition of the stochastic matrix p, it appears that each step (from state a to 
h) can be decomposed into two stages. The first one, corresponding to the stochastic matrix 
q, consists in the random choice of a trial move. Thus, q{Ca — > Ch) stands for the probability 
to perform a trial move from state a to h. Note that the Metropolis method requires q to be 
a symmetric matrix, so that the reverse move from state 6 to a has the same probability as 
the forward move. The next stage is the decision to either accept or reject this trial move, 
the probability of accepting or rejecting a trial move depending on the given weight function 
W{C), as shown in (plQl). If the trial move is rejected, the system stays in the same state, that 
is 6 = a. 

Now, forming a long enough random walk following this procedure, we have a new Monte 
Carlo estimator for the state density 

^[/.J..)c.gf '"^-^{^^^■^- . (12) 

where the Ci are the random configurations along the random walk of length A^. Here, the 
scale factor is defined by 

S={^W{C)^IN. (13) 

The second factor in the right-hand side of equation ([T2|) represents the (corrected) state density 
in the Monte Carlo sample, i.e., C(Jsmp(f^, J-, tt). Thus, each sampled configuration Ci is given an 
importance l/W{Ci) that is inversely proportional to its sampling probability P{Ci). 



2.3 Implementation of the Monte Carlo calculation 

In our calculation, each state belonging to the random walk corresponds to a configuration 
C, i.e., it is associated with a set of occupation numbers {n'^{C)]nl{C)} for all the neutron 
and proton single-particle states. Thus, the successive configurations are chosen randomly 
(according to the stochastic matrix p) as follows. First, one selects a nucleon at random 
(either a neutron with probability N/A or a proton with probability Z/A) and moves it to 
a randomly chosen available single-particle state. This generates a new configuration. (This 
single-particle trial move corresponds to our definition for the underlying matrix of the random 
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walk q, and is the most common choice for the elementary move. This sampling procedure is 
shown in ref. [^] to be generally more efficient than collective moves, with several particles 
moved simultaneously.) Then, one accepts or rejects this single-particle trial move according 
to equation (|T0|). This acceptance criteria will thus depend on W{C), which will be in general 
a function of Uc, Jc, and tiq- In our calculation, an energy dependence has proven to be 
sufficient, so that we choose W{C) as a function of Uc only. This weight function is arbitrary 
in principle: different choices give the same result for an infinitely long random walk. However, 
a judicious choice makes the method much more efficient than the simplest choice W{C) = 1. 
Indeed, we have to compensate for the rapid increase of the number of states with energy by 
selecting a correspondingly decreasing weight function. An appropriate choice would be to take 
W{C) proportional to the inverse of a Bethe-type law, 

W{C) = {a Ucf^^ exp [-2 (a Ucf'^] , (14) 

in which the parameter a is adjusted to ensure a more or less uniform distribution of the sampled 
states in the considered energy interval. We took this definition of W{C) in our previous work 
(see [^), for simplicity. In this paper, we prefer another definition based on a recursive model 
calculation. This will be discussed later on. 

Let us consider now the energy, spin, and parity assignment for each randomly sampled 
configuration. The excitation energy Uc of each configuration C (when assuming a system 
of non-interacting fermions) is simply given by the sum of the energies of all the occupied 
single-particle states 

Uc = Y.^l- + E 4 ■ nl{C) - Eo , (15) 

k k 

where the n'^'^{C) are the occupation numbers for the configuration C, and the e'^'^ are the 
single-particle energies for neutrons and protons. The ground state energy is noted as Eq. 
It is clear that, with our procedure, each excited level will be sampled with a probability 
proportional to its statistical weight, as expected. Suppose for instance that we are interested 
in the level corresponding to a (ld5/2)^ configuration, having a degeneracy d = =20. It will 
obviously be sampled with a probability proportional to the number of possible arrangements 
of 3 particles into 6 orbits, that is to d. The possible inclusion of pairing energy using BCS 
theory (removing this degeneracy) will be examined below. 

The spin Jc of each configuration C is chosen randomly among all the possible spins of 
its associated excited level; we choose the sampling probabilities of all the spins according to 
their respective degeneracies. In order to achieve this random choice, we proceed as follows. 
First, we determine the total spin of identical nucleons in each subshell. In order to account for 



the Pauli principle, we use the so-called Mayer- Jensen table [^, ^ which yields the number 
of times each spin occurs. Thus, the spin is chosen at random with a probability distribution 
determined according to this table, as explained in Appendix Then, the spins of the different 
subshells are coupled randomly, for either type of nucleons separately. For this purpose, we use 
the relation 

■J = ji+ h - min(rfi, rfa) (16) 
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for the coupling of two spins (ji and ^2) to a total spin J, where di and d2 are (discrete) 
random variables uniformly distributed in [0,2ji] and [0,2^2] respectively. It can be shown (see 
Appendix |^) that this relation yields a total spin J obeying the usual coupling rule of angular 
momentum for independent particles, with adequate weights (i.e., the sampling probability of 
each spin is proportional to its degeneracy). Thus, relation (p!6[) is applied recursively in order 
to couple the spins of all the subshells. Finally, the total neutron and proton spins are coupled 



at random, also using equation ([T6|) , to give the resulting spin of the configuration Jq- Note 
that, when doing a direct counting of the levels, the whole ensemble of possible spins (for each 
excited level) has to be enumerated. On the contrary, the Monte Carlo determination of the 
spin distribution is far simpler: one has just to choose one of these spins with the appropriate 
probability, which is much less time-consuming. 

The parity of the sampled configurations is also easily obtained by combining the parities 
of all the occupied single-particle states. Thus, the resulting parity ttc of a configuration C is 
given by 

= n K]"^^^^ n Kr^'^^ , (17) 

k k 

where the tt^''^ are the single-particle state parities. 

Finally, one has to determine the scale factor S needed to obtain an absolute state density. 
To this end, we calculate the total number of states A^rec up to our maximum energy ?7max (the 
upper limit of the energy interval of interest) by means of the Williams recursive method ||25|| : 

iVrec = / UJ,UU) dU . (18) 







Thus, we have a constraint on the Monte Carlo cumulative state density at energy 



max) 



A^(f/max) = i EE ^(U, J, vr) dU = K,, , (19) 







which gives obviously for the scale factor S = Nj-^c/N. This simple recursive method can neither 
yield the spin-parity distribution, nor take into account a residual interaction. However, we 
just need it to normalize our Monte Carlo sample state density UsmpiU, J, vr). Moreover, as it 
provides the whole state density curve as a function of the energy, Uj-edU), we can use it in 
order to define our weight function as explained in Appendix ^ This choice for W{C) ensures 
a more or less uniform distribution of the sampled states in the considered energy interval. As a 
consequence, the resulting state density uj{U, J, n) is more or less equally accurate in the whole 
energy interval. Note that we reject the sampling of excited states above the upper limit t/max 
by putting W{C) = for Uc > f/max- This recursive calculation of the state density takes a 
very short computation time. 

For illustrative purposes, we have calculated the Monte Carlo level density of ^^Fe up to 
an excitation energy of 30 MeV, using the realistic set of single-particle levels derived from 
a spherical Hartree-Fock calculation based on a Skyrme interaction In order to obtain 



a reasonable accuracy, the size of the sample A^ (which crucially determines the computation 
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time) is adopted equal to 10^. Figure [1| shows the derived cumulative state densities (i.e. the 
number of states up to a given energy U). The solid line represents the total cumulative state 
density (with both parities), while the positive- and negative-parity cumulative state densities 
are represented by dashed and dash-dotted lines, respectively. The curves show an overall 
behaviour in agreement with the statistical results, i.e. an asymptotic parabolic shape versus 
energy (in logarithmic scale). However, superimposed small-scale oscillations are pronounced 
at low energies, and disappear only at higher excitation. We also notice that the classical 
assumption of equipartition of both parities is only valid above ~ 15 MeV. Note that we 
obtain a high degeneracy for the ground state because of the neglect of residual interactions. 
Of course, this degeneracy would be removed if residual interactions were included in our model 
(see Section ^). For information, we also plot in Figure |I| the cumulative number of sampled 
states along the Monte Carlo simulation. It is evident that the excited states have been sampled 
approximately uniformly over the whole energy interval [0-30 MeV], as expected, so that the 
resulting accuracy is roughly constant (see Appendix 



3 Results 

We will show in this Section that the Monte Carlo procedure is very efficient for calculating 
an exact - at least in principle - level density with spin and parity dependence, whatever the 
mass number of the nucleus or the considered excitation energy. Other quantities, such as 
the fixed particle-hole number level densities, could also be obtained by this method. In this 
respect, the Monte Carlo method must be considered as a way to implement the calculation of 
a combinatorial level density, replacing a direct counting procedure when it is unfeasible (for 
large shell model spaces). However, a direct counting of the levels can be helpful when the 
number of excited levels is directly computable, in which case a Monte Carlo sampling would 
introduce undesirable statistical errors. 

In ref. [p^, we reported on the Monte Carlo derived level density of ^^Fe, ^"^Pb, and ^^°La. 



We chose the same size for the sample, i.e., N = 10 , which gave essentially the same accuracy 
in the three cases although the fraction of the whole set of excited states that is actually sampled 
is about 3 x 10"^ for ^'^Fe, 10"^ for ^ospb, and only about 2 x 10"^ for i^°La. With such a 
choice, it is remarkable that the Monte Carlo simulation takes about the same computation 
time for the three nuclei, although their mass numbers are very different. Some other cases can 
be found in ref. [0. In the following, we will limit ourselves to check the method by comparing 
with other traditional methods for calculating level densities. New results will be discussed in 



greater details in ref. ||30 |. 



3.1 Comparison with a direct counting 

In order to verify the reliability of our Monte Carlo procedure, we have performed calculations of 
the level density and spin-parity distribution via a direct counting of the levels (see ref. [0). We 
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have considered in particular the case of ^^Fe up to 30 MeV, neglecting the pairing interaction 
and using the single-particle levels from ref. |4^. Note that a combinatorial calculation at 
such excitation energies is still feasible in this mass number region, although at the cost of 
an impressive cpu time. The total state density derived with the Monte Carlo method could 
hardly be distinguished on Figure |l| from the direct counting , so that we did not plot the 
latter. In fact, the tiny difference is simply due to the statistical noise (which can be made 
arbitrarily small by increasing N). Of course, the recursive method also gives identical results, 
but, nevertheless, it is impracticable when considering residual interactions. 

Let us turn now to the spin dependence of the ^^Fe level density. In Figure |^, we plot the 
spin distribution of the levels in the energy interval [20-21 MeV] derived by the Monte Carlo 
method, along with the direct counting calculation . The square symbols and the solid line 
represent, respectively, the result of the Monte Carlo method and the direct counting. The 
Monte Carlo spin distribution is obtained by putting W{C) = 1 from f/ = to the upper limit 
U = 21 MeV. Thereby, the states with high excitation energy are sampled preferentially, and 
the accuracy on the spin distribution is increased. The adopted sample size is = 50 x 10^, 
in order to get very good statistics (to have negligible statistical errors). Of course, such an 
accuracy might be unnecessary for practical applications, but our intention here is just to 
show that the Monte Carlo result is asymptotically exact. Indeed, both Monte Carlo and 
direct counting calculations are in very good agreement, showing the validity of the random 
spin coupling procedure explained in Section ^ Note that, in this high energy domain, the 
number of excited levels is as high as 2 x 10^, so that the spin distribution is expected to be 
consistent with the asymptotic distribution of the statistical model (see Figure 0). The positive- 
and negative-parity spin distributions obtained by the Monte Carlo method are also shown in 
Figure 0. It appears that, in spite of this large number of levels, the parity equipartition is 
not achieved yet in the considered energy interval, illustrating the need for an appropriate 
treatment of parity (see [^ ). 

In the conventional statistical model, the dependence of the state density upon the angular 
momentum projection M has the gaussian form R ffl 



u;{U,M) = u;{U) 



exp(-MV2a2) 
V27ra2 



(20) 



where cx^ is the spin cut-off parameter, and u!{U) stands for the total state density. The resulting 
spin-dependent level density can be written as P| 



p{U,J) = uj{U,M = J) -uj{U,M = J + 1) 
d 



dM 
2J + 1 



uj{U, M) 



2(271)1/2^3 



exp 



. M=J+l/2 

(J +1/2)^ 
2a2 



uj{U) 



(21) 



where J denotes the total angular momentum of the levels at excitation energy U . 
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We check in Figure ^that our Monte Carlo derived spin distribution for ^^Fe in [8-10 MeV] 
looks similar to the gaussian-shaped asymptotic distribution. We plot the M-distribution of the 
excited states (i.e., u{U,M) as a function of M). The agreement with the gaussian form (|2y) 
is evident, as expected given the excitation energy. We estimate an effective spin cut-off pa- 
rameter 0"^ by calculating the variance of this M-distribution. This yields the value cr^ ~ 18.8, 
corresponding to a maximum in the J-distribution at a spin around a ~ 4. The gaussian 
distribution with this adopted value for cj^ is plotted. Figure || also shows the correspond- 
ing J-distribution of the excited levels simulated by Monte Carlo, as well as the asymptotic 
distribution ( pT]) obtained with the same value of a^. It proves that one can be confident in 
the asymptotic result if one is concerned with an energy interval containing a sufficiently large 
number of levels. Indeed, the number of levels amounts to 8 x 10^ in our case. However, when 
considering an interval containing fewer levels, the assumption of a gaussian distribution be- 
comes very poor (see e.g. |3^, |^). It is interesting to note that the discrepancy between the 



curves is more pronounced when looking at the J-distribution, showing that the latter is more 
sensitive to deviations from asymptotic results. 

Finally, the evolution of the parity distribution of the excited levels in ^^Fe as a function 
of energy (up to 30 MeV) is shown in Figure ^. We plot, as a function of U, the cumulative 
parity asymmetry defined by 

A = (22) 

N+{U)+N-{U) ' ^ ' 

with N^{U) = uj^{U) dU. The dotted line corresponds to the Monte Carlo results, whereas 
the solid line is calculated by a direct counting of the levels [Q. Once more, both curves are 



almost indistinguishable. This proves that the Monte Carlo method is a powerful alternative 
to the traditional combinatorial method; it reproduces the parity-asymmetry curve derived by 
direct counting within an asymptotically small statistical error and with a huge gain of cpu 
time. Note that the parity equipartition, generally assumed to be achieved at low energies, is 



only reached above ~ 15MeV. This point is examined in greater detail in ref. |42 



3.2 Comparison with a microscopic statistical model 

In order to illustrate the reliability of our Monte Carlo procedure, we have also compared 
our results with the predictions of a numerical statistical model. This model, based on the 
traditional methods of statistical mechanics, takes into account the shell effects by using a 
realistic discrete single-particle spectrum (see e.g. 0). The total energy E, the number of 
neutrons N and protons Z, and the entropy S are determined from the single-particle energies 
el (for neutrons) and (for protons) by the relations 
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S = n- onN - azZ + (3E , (26) 

with 

Q = \nZ = Y. ln[l + exp(a^ - /3e^)] + ln[l + exp(az - f^el)] , (27) 

k k 

where Z is the partition function, and oat, a^, [3 = 1/t are the Lagrange muhiphers related, 
respectively, to the number of neutrons, the number of protons, and the total energy. One also 
defines the quantity D as the 3x3 determinant of the second derivatives of VL with respect to 
Oat, ttz, and /3, evaluated at the saddle point. All these quantities can be calculated numerically, 
and the resulting state density for a two-component system is given by 

^(^>^>^)= (2^)?2^i/2 - (28) 



These calculations were performed using the statistical method from ref. |^9[, based on a realistic 
spectrum of single-particle levels from ref. [^], neglecting the pairing interaction. 



In Figure |^a, we show the state density for ^^Fe calculated by this statistical method for 
comparison with the combinatorial (Monte Carlo) result. The solid line represents the state 
density obtained with our Monte Carlo method, plotted in bins of 100 keV, while the dashed line 
corresponds to the statistical model calculation. Both curves were calculated using the same set 
of single-particle levels. Note that the noisy behaviour of the Monte Carlo state density at low 
energies comes from the fact that a lot of bins are devoid of levels. Of course, the curve could be 
smoothed by increasing the bin width, but the oscillations (due to the discrete structure of the 
single-particle spectrum) would then be smeared out. The agreement between statistical and 
combinatorial state densities is reasonable, even if those oscillations disappear completely in 
the statistical formalism. The statistical state density can thus be viewed as an average of the 
combinatorial state density, which is obviously a consequence of the grand-canonical approach 
(the energy is only fixed on average). Note that the statistical curve looks higher than the 
averaged combinatorial curve, but this is simply due to the logarithmic scale. It has to be 
stressed that the discrepancies between both methods (ranging up to one order of magnitude) 
might be significant e.g. when calculating reaction cross sections with statistical or Monte 
Carlo level densities. We also consider the case of a heavier odd-odd nucleus in Figure We 
compare the combinatorial (Monte Carlo) and statistical state densities for ^^°La (assumed to 
be approximately spherical, and without pairing). The solid and the dashed line correspond, 
respectively, the Monte Carlo and the statistical state density. The latter curve also appears as 
an averaged combinatorial state density, even if a little systematic error seems to be present. 
A small statistical noise is still visible in the Monte Carlo state density, related to the choice 
of the sample size N = 10^. 

Let us consider now the spin distribution in this microscopic statistical model, and in par- 
ticular the evolution of the spin cut-off parameter with energy. Within the framework of the 
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statistical model, one can generalize the formalism by adding a supplementary Lagrange multi- 
plier 7 related to the angular momentum projection M (see [^). It can be shown that, except 
for extremely large values of M, the Lagrange multiplier 7 is sufficiently small to obtain a 
good approximation by expanding the statistical expressions in powers of 7, and by keeping 
the terms up to 7^ only. The thermodynamic quantities of interest can be written in terms of 7 
as in equations By eliminating 7, one gets an expression for the M-dependent state 



density similar to equation (^Of) 



iV, Z, M) = AT, ^) ^^P( MV2a^) ^^9) 

where uj{E, N, Z) is the total state density defined in equation (|28|). The spin cut-off parameter 
which determines the width of the M-distribution is given by the expression: 

^' = \ E«)'sech2^(/3e^ + \ E(K)'sech2^(/3e^ - az) , (30) 

k k 

where and are the single-particle magnetic quantum numbers for neutrons and protons, 
respectively. The resulting spin-dependent level density p{E, N, Z, J) can be expressed as in 
equation (pT]). 

In Figure |^, we compare the effective spin cut-off parameters of three nuclei (^^°Sn, ^^^Dy, 
and ^''^Pb) derived from our combinatorial (Monte Carlo) method and from a numerical statis- 
tical model (ref. [^, |5^). Both predictions make use of the same Woods-Saxon single-particle 
level scheme. The Monte Carlo spin cut-off is derived by calculating the variance of the M- 
distribution, as previously. As can be seen, the Monte Carlo method leads to a rather good 
prediction of the energy- dependent effective spin cut-off parameter, showing that one can be 
confident in the method. Note that, even if this agreement for is very satisfying, the sta- 
tistical assumption of a gaussian distribution is not always justified (see 0), so that our 
Monte Carlo procedure can become essential in some cases. 



4 Pairing interaction 



In our previous work [0, we restricted ourselves to the Monte Carlo evaluation of the level 
density without the inclusion of pairing force. In this Section, we explore the possibility to ac- 
count for the pairing interaction in the Monte Carlo method. The pairing interaction is usually 
taken into account in combinatorial calculations by use of the BCS theory (see ref. p2| , ^). 
For example, Hillman and Grover corrected their combinatorially calculated level density 
for pairing forces by solving the BCS equations (at zero temperature) for each excited config- 
uration. They used two so-called "interacting odometers" to scan the nuclear configurations 
systematically, one for the unpaired nucleon configurations (determining the spin distribution), 
the other one for the pair configurations. They extended the "blocking" method developed in 



ref. ||51|] for treating odd-A nuclei (with one unpaired nucleon) to the case of configurations 
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containing more than one unpaired nucleon: all orbital pairs containing unpaired nucleons are 
blocked (i.e., made unavailable for pair diffusion by being omitted from the BCS calculation). 
In their model, the total energy of a (proton or neutron) configuration C is given by 

^c = E'e^ + E"2e.<c-Ac/G', (31) 

k k 

where G is the pairing strength parameter, and v\ q is calculated from 

""^'^ [(^.-Ac)^+%]^/4 ' ^^^^ 

Ac and Ac being obtained by solving the pair of BCS equations 

E"2<c = ^c, (33) 

k 

E"[(e.-Ac)^ + A^]-V2 = 2/G. (34) 

k 

Here, Ac and Ac are the gap parameter and Fermi energy, respectively, while r/c stands for 
the number of paired nucleons. The summation is made over orbitals containing unpaired 
nucleons, while Y^' is made over orbital pairs in which there is not an unpaired nucleon. Note 
that, as the excitation energy increases, more and more orbital pairs are blocked, so that only 
the trivial solution A = can exist in general above a certain energy. This behaviour is in 



qualitative agreement with the result of the temperature-dependent BCS theory |T^ . 

However, the model of Hillman and Grover has an inherent drawback related to the con- 
figurations involving promoted pairs. Due to its variational nature, the BCS method cannot 
treat a general configuration of pairs, but only yields the energy of the particular configuration 
where the pairs are all placed in the lowest available (unblocked) orbitals. They accounted 
for this problem by adding to the BCS energy the difference in single-particle energy of that 
particular configuration of pairs and the one under consideration. They claimed that, even if 
this is clearly a rough approximation, the resulting error is not very serious because the fraction 
of levels with pair excitation is small. 

The combinatorial model of Herman and Reffo is slightly different, although the pairing 
interaction is also included by applying the BCS theory to each configuration. As they were 
interested in estimating a level density with fixed exciton numbers, they simply neglected the 
excitations involving promoted pairs. All the excited particles are taken as noninteracting 
excitons, and the orbitals in which they are placed are excluded from the BCS consideration, 
even if two particles occupy time-reversed orbitals. They also justified this approximation by 
arguing that the number of such configurations involving promoted pairs is very small. 

In conclusion, it is our opinion that this problem of promoted pairs is due to the hybrid 
character of those methods. Indeed, BCS is inherently a statistical theory, and is therefore 
incompatible with any combinatorial procedure for estimating the level density. In fact, only 
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a kind of "combinatorial" treatment of the pairing interaction would give a consistent solution 



to the above problem. This is our purpose in ref. ||52[, where we propose a quantum Monte 
Carlo method to treat exactly the pairing force in nuclei. This method in its present status is, 
however, limited to the properties of the ground state, and is not directly applicable to the level 
density problem. Thus, in this paper, we have to limit ourselves to an approximate inclusion 



of pairing by use of the standard BCS theory, close to what is done in ref. |2^, and we show it 
is feasible at essentially no cpu cost. 

For each sampled excited configuration, we solve the BCS equations in order to obtain 
its pairing energy Pc- We also use the so-called blocking method that is the orbits 

occupied by unpaired nucleons are blocked (i.e., unavailable for pair diffusion). However, we 
treat the problem of the promoted pairs in a slightly different way: first, the number of unpaired 
nucleons occupying the Fermi level is taken equal to the seniority s derived from our random 
spin coupling algorithm, as explained in Appendix |B[ This will ensure that the spin attribution 
of the low-lying levels is coherent with the single-particle level scheme. Second, all the excited 
particles (and created holes) on the other levels are taken as noninteracting excitons, even if 



two particles (or holes) occupy time-reversed orbitals, in analogy with what is done in ref. ||24 
We think that this is the simplest approximate prescription to properly construct the first 
excited states, while neglecting excitations involving promoted pairs at higher energies. Thus, 
the energy Uc of each sampled configuration is corrected for pairing and replaced by Uc — Pc 
in equation (|12|), when computing the Monte Carlo state density: 



up{U,J,7i)c^S}_^ ^^j^^j^^ , (35) 

where Pi stands for the pairing energy of the configuration Ci (that is the sum of the pairing 
energies for neutrons and protons). Note that the energy used to calculate the weight function 
(see Appendix 0) does not include this pairing term, so that the state density without pairing 
can also be computed at the same time. This state density is useful since it allows to perform the 
normalization according to c^rec, which do not incorporate pairing energy, as already mentioned. 
In fact, the Monte Carlo simulation proceeds as if there was no pairing (the weight function is 
inversely proportional to uj^-ec), and then the actual sampling of the configurations corrects for 
the effect of pairing. 

In order to illustrate the effect of pairing, we have considered the isobars ^^^Nd (even-even) 
and i^^Pj, (odd-odd). For both nuclei, we use the spherical single-particle level scheme resulting 
from a Hartree-Fock + BCS calculation from ref. [J^. We solve the BCS equations for all the 



excited configurations with the values for the pairing strength parameters Gn = 2.25/A^°'^ MeV 
(for neutrons) and Gp = 2.00/Z°-^ MeV (for protons) derived in ref. |^6| for the ground state. 



In Figure ^ we plot the cumulative total state density including pairing effects derived with our 
Monte Carlo method for both ^^^Nd (solid line) and ^^^pj. (dash-dotted line), as a function of 
the excitation energy U. The effect of deformation is neglected here. The dashed and the dotted 
lines represent the cumulative state density for both nuclei in the absence of pairing. The latter 
curves clearly illustrate the existence of shell effects: as ^^^Nd is situated at a neutron shell 
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closure, its level density is lower than that of the neighbour nucleus ^^^Pr. In fact, given the 
similarity between the corresponding single-particle level schemes, this can be simply described 



as a Rosenzweig effect |]S3[ , resulting in a shell-correction shift in the energy scale. In addition, 
one can see in Figure |^ that, for both nuclei, the pairing force is responsible for a shift of 
the level density curve to the right, according to the classical phenomenological models (see 
e.g. [0). Roughly speaking, this shift corresponds to the condensation energy of the ground 
state due to pairing. Therefore it is clear that, due to the even-even character of ^^^Nd, its 
level density has to be shifted more than that of ^^^Pr. (The condensation energy is found to 
be about 4.6 MeV for ^^^Nd, and 3.3 MeV for ^^^Pr.) The difference between both shifts is 
approximately equal to A„ -|- Ap, with A„ and Ap being respectively the gap parameter for 
neutrons and protons in ^^^Nd. The Monte Carlo method is, of course, more accurate than a 
simple phenomenological model since it yields a shift varying with energy, and thus provides a 
supposedly more realistic level density. 

Note that the extra amount of cpu time needed for the inclusion of pairing in our Monte 
Carlo procedure was found to be only about 30-40%. This fact seems surprising at first glance, 
since it is well known that a combinatorial computation becomes in general extremely time- 
consuming when treating residual interactions. It originates from the fact that only about 
1% of the Metropolis trial moves are accepted on average along our random walk, so that the 
extra cpu time needed to compute the pairing energy (by solving BCS equations) only concerns 
about 1% of the configurations. Accordingly, the energy of the remaining (99%) configurations 
does not need to be calculated, allowing a significant gain of computing time. The latter 
configurations are necessary, however, to associate the appropriate statistical weights to the 
former 1% configurations that are effectively sampled. 

Therefore, we think that a Monte Carlo evaluation of the level density with inclusion of 
pairing is very promising. If the pairing force could be accounted for in a consistent way (i.e., 
without the use of BCS theory), it would be the only available method for solving exactly this 
problem. It would not depend on the grand-canonical statistical formalism (finite-temperature 
BCS equations [|1^]) and its associated difficulties, such as the spurious phenomenon of a sharp 



superfiuid- normal phase transition in nuclei |54]. The exact solution of this problem via a 



Monte Carlo procedure is the subject of a future work. 



5 Conclusion 

The Monte Carlo method is shown to be a reliable technique for calculating nuclear level den- 
sities with a sufficient accuracy within computation times that are much shorter than those 
demanded by traditional combinatorial methods (i.e., a direct counting). Also, owing to its 
rapidity, this method makes possible the estimate of the shell-model uncertainties (i.e., those 
associated with the single-particle level scheme). In the absence of residual interactions, how- 
ever, the recursive method is exact as well and requires approximately the same amount of 
computation time. In any case, it is necessary to calculate the Monte Carlo scale factor by the 



19 



recursive method. Thus the Monte Carlo method really becomes a powerful tool when study- 
ing nuclear level densities with pairing interactions, as explored in Section ^. It would even be 
possible to include neutron-proton interactions in this method, which cannot be considered by 
any other method. 

The Monte Carlo approach also turns out to be worthwhile to predict the excited levels 
spin and parity distributions. Considering the crude approximations made to obtain those 
distributions in the framework of the statistical model, the proposed method is an interesting 
alternative to the usual statistical methods used in nuclear reaction calculations. Moreover, 
the technique outhned in this paper could be extended to the case of deformed nuclei. More 
detailed calculations and a confrontation with experimental data are presented in a subsequent 
article 
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Appendices 



A Statistical errors 

Due to the principle of the Monte Carlo method (i.e., the sampling of a finite number of 
configurations), it is clear that a statistical noise is superimposed on the derived total state 
density and the spin-parity distribution. We will only consider here the case of the total state 
density, but the same reasoning also applies to the calculation of the statistical error on other 
quantities. 

The continued evolution of the random walk yields successive configurations C that are 
distributed according to the weight function W{C). This distribution is then sampled by 
repeated application of the stochastic process p{Ci — > Q+i) defined in equation (0). The 
resulting Monte Carlo estimator for the total state density is thus given by 



u{U) 




P(C) 

(36) 



where S = J2c W{C)/N is the scale factor, and E[ ] means the expectation value for a random 
sampling of the configurations C according to the probability distribution P{C) oc W{C). 
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As explained in Section |2.2| , the idea of the importance samphng procedure is to achieve a 
more or less uniform sampling in the considered energy range. Let us consider an energy bin of 
width Af/, situated at an excitation energy U, containing a non-zero number of states. (If the 
bin is devoid of states, it must simply be excluded from the following considerations.) Let us call 
{C} the ensemble of configurations belonging to this bin (i.e., for which U < Uq < U + AU), 
and let us give them the same weight W{U) for simplicity. Assuming that the state density 
is constant inside this bin and equal to Cd{U), the weight of that bin is thus approximately 
W{U)u{U)AU. The condition of a uniform sampling imposes thus that 

W{UMU)AU/ J2 ^iC) ^ AU/U^^, , (37) 
c 

where f/max is the upper limit of the considered energy interval. In practice, since uj{U) is 
unknown, it is replaced by the derived state density from the recursive method. Equation (|37|) 
corresponds to our definition of the weight function W{C). 

According to equ. (^6]), the estimator for the number of states in that bin is obviously 



N 



{S/W{U)) Yl - U)Q{U + AU-UcJ, (38) 



i=l 



where G(x) stands for the Heaviside step function. By integrating (|36| ) between U and U + AU, 
the expectation value of MiJJ) is trivially given by 

E[7V([/)] = / u{U)dU = uo{U)AU . (39) 

In order to estimate the statistical noise on J\f{U), let us define first the probability of choosing 
a random configuration in this bin as "P = W{U)uj{U)AU/ J2c^{C) ~ AfZ/f/max- Assuming 
that the configurations are independent (see below), one has for its variance 

Var[7V(t/)] = [S/W{U)y NV{1 - V) . (40) 
Thus, the (relative) statistical error on N'{U) is simply given by 

,/^- (41) 

If the sampled state distribution is more or less uniform, V will be approximately equal to 1/-B, 
B being the number of bins in the energy interval. Thus, the (relative) statistical error will be 



'l-V 



on the order of \JB/N, i.e. the inverse of the square root of the number of states which have 
been sampled in that bin (whatever the actual number of states in that bin). 
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Finally, another fact has to be taken into account when estimating the statistical error on 
various quantities. Indeed, the successive configurations along the random walk are clearly not 
statistically independent since each one is generated by moving one particle at most from the 
previous one; that is, Cj+i is likely to be in the neighbourhood of Ci even if the configurations 
are distributed properly as the walk becomes long. As a result, the above expression of the 
variance is clearly underestimated. This can be quantified by calculating the auto-correlation 
function of some estimator, e.g. the energy Uq- In practice, the auto-correlation length Lcorr can 
be computed as the number of random steps for which this function becomes reasonably small. 
Then, when estimating the statistical errors on any quantity, the number of sampled states 
has to be replaced by N/Lcon, measuring the number of independent sampled configurations 
during the Monte Carlo run. 



B Implementation of the random spin coupling 

The calculation of the spin distribution is simplified since we restrict the method to spherical 
nuclei here. The orbital angular momentum j becomes thus a good quantum number, allowing 
us to define subshells. The number of levels having a given total angular momentum J for a 
particular configuration of subshells can be calculated by finding the number of ways the com- 
ponents j's can couple to J according to usual coupling rules. We show in i) how we achieve 
this coupling in a random manner. However, when we treat two (or more) identical nucle- 
ons in a same subshell, the restrictions imposed by the Pauli exclusion principle for fermions 
complicate this coupling, as explained in ii). Note that we do not have to construct explicitly 
the eigenf unctions of by use of the Clebsch-Gordan or fractional parentage coefficients, as 
the method only serves to determine the correct number of states with their good spin. As 
a matter of fact, we make an arbitrary correspondence between the \m,im,2 ■ ■ ■) states and the 
\JM) states. We simply choose at random the |mim2 ■ ■ ■) states so that the corresponding 
\JM) states are sampled with the correct statistical weights. Thus, we do not calculate an 
intermediate M-dependent state density uj^U, M) in order to deduce p{U, J) through the usual 



relation p(f/, J) = uj{U^ M = J) —u}{U, M = J + 1), as in ref. p3|. Rather, we calculate directly 



the spin distribution, which allows to properly incorporate pairing effects. 

i) Spin coupling for different orhitals: 

Let us show how to couple randomly the spins of the different subshells (for neutrons or 
protons). As mentioned in Section |2]^, we use the following relation for the coupling of two 
spins (ji and j2) to a total spin J: 

J = ii+ J2 - min((ii, ^2) , (42) 

where di and d2 are discrete random variables uniformly distributed in [0, 2ji] and [0,2^2], 
respectively. This relation yields a total spin J obeying the usual coupling rule of angular 
momentum for independent particles, with adequate weights (i.e., the sampling probabihty of 
each spin is proportional to its degeneracy). 
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In order to prove this, let us calculate the probability distribution of the random variable 
J defined by equation (^21) . Since di and c?2 are uniformly distributed, they have the following 
cumulative distribution functions: 

F,{d,) = IItTT if0<^i<2ji, (43) 
I 1 if rfi > 2ji , 

U if Ci2 > 2j2 . 

Let us define the random variable d = min((ii, ^2), and let us determine its cumulative distri- 
bution function F{d) using 

F{d) = Prob [min(c?i,rf2) < c?] 

= Prob [{di < d) or (c/2 < d)] 

= 1 - Prob[di > d\ ■ Prob[rf2 > d\ 

= F,{d) + F2{d) - F,{d) ■ F2{d) , (45) 

where we have used the fact that di and d2 are independent random variables. Thus, if d 
reaches the maximum value of either di or d2 (i.e., 2ji or 2^2 respectively), one has Fi{d) = 1 
or F2{d) = 1, so that 

F{d) = 1 a d> min(2ji,2j2) • (46) 

This implies that the probability of having d > min(2ji, 2^2) is equal to zero, so that J cannot 
be less than ji + j2 — min(2ji, 2j2) = \ji — ^2!, as expected. Using equations (^3[), (^) and (|i5|) 
we can calculate the probability distribution f{d) of the random variable d as 



f{d) = F{d)-F{d~l) 
2{3i+32-d) + l 



(2ji + l)(2j2 + l) 
Thus, the probability distribution of J is simply given by 



with < d < min(2ji, 2j2) • (47) 



HJ) = /(J1+J2-J) 
(2J+1) 



(2ji + l)(2j2 + l) 



with IJ1-J2I <J<Ji+j2 , (48) 



which exactly corresponds (within a multiplicative factor) to the statistical weight of the spin 
J. Note that the denominator comes from the normalization condition 

ii+i2 

E <^>(^) = l (49) 

J=\jl~j2\ 
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of the probability distribution $(</). This procedure can be obviously generalized for the 
coupling of n independent spins by applying relation (^21) recursively. 



ii) Identical nucleons in a single orbital: 

For identical nucleons in a subshell, only those states are permitted for which the Pauli 
exclusion principle is obeyed. Thus, the Monte Carlo implementation of the spin coupling 
is not straightforward, as for different orbitals (case i). A rapid method for handling this 
problem is to put the required spin distribution in a table. This table, constructed according 
to the Mayer- Jensen table |^, yields the number of times each spin J occurs for various 
configurations (j)'^, where j is the spin of the subshell and k is the number of identical nucleons 
occupying the subshell. 

As an example, we present in Table |l| a list of the possible angular momenta J for an orbital 
j = 9/2. Note that the multiplicity of the states due to the degeneracy for M is taken into 
account. The states are labelled by the seniority quantum number s (corresponding to the 
number of unpaired particles in the orbital), with 

s = 0, 2, 4, • • • min(A;, 2j + 1 — k) for k even , 

s = 1, 3, 5, ■ ■ ■ min(A;, 2j + 1 - k) for k odd . (50) 

For instance, the various spins of the two-particles states {k = 2) are obtained by summing 
the rows s = and s = 2 in Table |I]. The last column in Table |l| shows the total multiplicity 
of the states of seniority s. Note that when the subshell is more than half filled (i.e., when 
k > j + 1/2), one applies the same procedure with the number of holes 2 j -|- 1 — k. Our Monte 
Carlo procedure to give the spin of a configuration (j = 9/2)^ consists in choosing at random 
one of the values J=0, 2, 4, 6, 8 with the corresponding probabilities 1/45, 5/45, 9/45, 13/45, 
17/45. The total weight (needed to normalize the probabilities) is simply 1 -|- 44 = 45 = (^2^ 
At the same time, the seniority s is chosen equal or 2 according to the probabilities 1/45, 
44/45. The value of the seniority s yields the number of unpaired particles in the considered 
subshell, needed when treating the pairing interaction (see Section |^). The same table for 
different values of j has been computed, and an analog procedure is used to yield a (random) 
spin for all the orbitals. 
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Table 1: List of the statistical weights of the possible states for different values of spin J 

and seniority s in a j = 9/2 subshell. The upper panel corresponds to states with an even 
number of particles A; (s = 0, 2, • • • A;), while the lower panel corresponds to states with odd k 
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Figure captions 



Figure 1: Cumulative state density for ^^Fe. The solid line represents the total (with both 
parities) density, while the positive- and negative-parity densities are represented by dashed 
and dash-dotted lines, respectively. The dotted line corresponds to the cumulative number of 
sampled states. The adopted size of the sample is iV = 10^. 



Figure 2: Spin distribution (normalized to 1) of the excited levels in the energy interval [20- 
21 MeV] for ^^Fe. The squares represent the results of our Monte Carlo procedure, while the 
solid line represents those of a direct counting calculation ESI. Both calculations were made 
using the same single-particle level scheme p6l. The triangles and the crosses correspond. 



respectively, to the Monte Carlo positive- and negative-parity spin distribution. A sample size 
AT = 50 X 10^ is adopted. 



Figure 3: Distribution of the excited state angular momentum projection M and of the 
excited level spin J for ^^Fe in the energy interval [8-10 MeV]. The squares correspond to 
uj{U, M) resulting from our Monte Carlo calculation, while the solid line represents the gaussian 
distribution from eq. (^). The triangles stand for the Monte Carlo derived p(f/, J), while the 
dashed line represents the asymptotic distribution from the statistical model (eq. ^l]). The 
derived value for the effective spin cut-off parameter is ^ 18.8. 
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Figure 4: Cumulative parity asymmetry A as a function of excitation energy U for ^^Fe. The 
dotted line represents the result of our Monte Carlo procedure, while the solid line represents 
those of a direct counting Note that, for the sake of emphasizing the small discrepancies 
between both curves, we have slightly shifted the dotted curve to the right. Both calculations 



were made using the same single-particle level scheme [46 



Figure 5: State density for ^^Fe (a) and ^"^^La (b). The solid line represents the combinatorial 
state density obtained with our Monte Carlo method, plotted in bins of 100 keV, while the 
dashed hne corresponds to the statistical model calculation using the same spectrum of 
single-particle levels (from [46|). The statistical curve appears as an averaged combinatorial 
state density. 



Figure 6: Comparison of the effective spin cut-off parameter cr^ as a function of energy U, 
for ^^°Sn, ^^^Dy, and ^°^Pb. The symbols represent the values derived from our Monte Carlo 
method (each symbol corresponds to a Monte Carlo simulation) , while the solid lines stand for 
the predictions of the statistical model Both evaluations have been obtained making use 



of the same Woods-Saxon single-particle level scheme. This figure is taken from 



Figure 7: Cumulative state density including pairing effects derived with our Monte Carlo 
method for two isobar {A = 142) nuclei, as a function of excitation energy U. The solid and 
the dash-dotted line correspond, respectively, to the case of ^^^Nd and ^''^Pr, while the dashed 
and the dotted lines represent the cumulative state density for both nuclei in the absence of 
pairing. We took the single-particle level schemes and the pairing strength parameters from 
ref. H (i.e., Gn ^ 0.10 MeV and Gp ^ 0.11 MeV for both nuclei). 
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